#ifndef EVALUATEMONOMIALS_H
#define EVALUATEMONOMIALS_H

#include <type_traits>
#include <cmath>
#include "Core/static_for.h"
#include "Core/numlib.h"

struct identityOfMonomials
{
  template <int Dim, class... Ts>
  static Real transform(const Vec<Real, Dim> &p, Ts... indices) {
    Vec<int, Dim> beta { decltype(indices)::value ... };
    Real r = 1.0;
    for(int d = 0; d < Dim; ++d)
      r *= std::pow(p[d], beta[d]);
    return r;
  }
};

struct gradOfMonomials
{
  template <int Dim, class... Ts>
  static Vec<Real,Dim> transform(Vec<Real, Dim> p, Ts... indices) {
    Vec<int, Dim> beta { decltype(indices)::value ... };
    Vec<Real,Dim> r;
    for(int d = 0; d < Dim; ++d) {
      if(beta[d] == 0) {
        r[d] = 0.0;
      } else {
        r[d] = 1.0;
        --beta[d];
        for(int k = 0; k < Dim; ++k)
          r[d] *= std::pow(p[k], beta[k]);
        ++beta[d];
      }
    }
    return r * beta;
  }
};

template <int Dim, int n>
class EvaluateMonomials
{
public:
  template <class T_DiffOp, class T_Mod, int D0, class... Ts>
  static void doEvaluate(const T_DiffOp &diffOp, const T_Mod &modifier, Real* &where2Fill,
                         const Vec<Real,D0> &p, const Vec<Real,D0> &nml, const Vec<Real,D0> &dx, Ts... moreArgs)
  {
    static_for<0, n+1>::execute([&, moreArgs...](auto _i_, int) {
      constexpr int i = decltype(_i_)::value;
      EvaluateMonomials<Dim-1, n-i>::doEvaluate(diffOp, modifier, where2Fill, p, nml, dx, _i_, moreArgs...);
      return 0;
    }, 0);
  }
};

template <int n>
class EvaluateMonomials<1, n>
{
public:
  template <class T_DiffOp, class T_Mod, int D0, class... Ts>
  static void doEvaluate(const T_DiffOp &diffOp, const T_Mod &modifier, Real* &where2Fill,
                         const Vec<Real,D0> &p, const Vec<Real,D0> &nml, const Vec<Real,D0> &dx, Ts... moreArgs)
  {
    static_for<0, n+1>::execute([&, moreArgs...](auto _i_, int) {
      *where2Fill += modifier(diffOp.transform(p, _i_, moreArgs...), nml) / identityOfMonomials::transform(dx, _i_, moreArgs...);
      ++where2Fill;
      return 0;
    }, 0);
  }
};

#endif //EVALUATEMONOMIALS_H
